function svl_jpe_tables()
% Section 5.2

clear
clc
close all

[ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters();

% table2(ell,theta,sigma,pi,rho,eta,i,M,A,B) %    5.346179 seconds
% table3(ell,theta,sigma,pi,rho,eta,i,M,A,B) % 3748.086483 seconds

%=========================================================================
% parameters
%-------------------------------------------------------------------------
function [ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters()
% parameter setting
ell   = 0.5;       % probability of being C-type
theta = 0.5;       % C-type's bargaining power
sigma = 1;         % CRRA IES
pi    = 0.95;      % 1 - default rate
rho   = 0;         % recovery rate
eta   = 0;         % IRS
i     = 0.1;
M     = 1;
A     = 0.05;
B     = 0.05;

%=========================================================================
% subfunctions
%-------------------------------------------------------------------------
function u = u(q,sigma)
% CRRA utility: u(q)
if sigma < 1
    u = q.^(1-sigma)./(1-sigma);
elseif sigma == 1
    u = log(q);
end

function mu = mu(q,sigma)
% marginal utility: u'(q)
mu = q.^(-sigma);

function w = w(q,sigma,theta)
% w function in effective bargaining power
w = theta + (1-theta)*mu(q,sigma);

%=========================================================================
% solvers
%-------------------------------------------------------------------------
function [G_eC,L_A,L_B,q0A,q1A,q0B,q1Bn,q1Bd,eNn,eNd]=partialEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given eC

lb = [0,0,0,0,0,0,0]; % lower bound constraint
ub = [1,1,1,1,1,1,1]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9,0.9,A/(A+B),A/(A+B)]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off'); 
if eC > 0 && eC < 1
    nonlcon  = @(x) fminconstr(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);
elseif eC == 0
    nonlcon0 = @(x) fminconstr0(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon0,opts);
elseif eC == 1
    nonlcon1 = @(x) fminconstr1(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon1,opts);
end

q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
q1Bd = x(5);
eNn  = x(6);
eNd  = x(7);

if eC > 0 && eC < 1
	a_CA_n = eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
	a_CB_n = (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
	a_CA_d = eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
	a_CB_d = (1-eNd)*(1-ell)/((1-eC)*ell+(1-eNd)*(1-ell))^(1-eta);
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	S_CB_d = theta*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
	L_A    = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
	L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1) + ell*(1-pi)*rho*a_CB_d*theta/w(q1Bd,sigma,theta)*(mu(q1Bd,sigma)-1);
	S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
	S_CB_tilde  = - i*q0B - L_B_rho_bar*((1-theta)*(u(q1Bn,sigma)-u(q0B,sigma))+theta*(q1Bn-q0B)) + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CB_d*S_CB_d);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / (pi+(1-pi)*rho);
elseif eC == 0
	a_CA_n = 0;
	a_CB_n = 1-ell;
	a_CA_d = 0;
	a_CB_d = 1-ell;
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	S_CB_d = theta*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
	L_A    = 0;
	L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1) + ell*(1-pi)*rho*a_CB_d*theta/w(q1Bd,sigma,theta)*(mu(q1Bd,sigma)-1);
	S_CA_tilde  = - i*q0A + ell*(u(q0A,sigma)-q0A);
	S_CB_tilde  = - i*q0B - L_B_rho_bar*((1-theta)*(u(q1Bn,sigma)-u(q0B,sigma))+theta*(q1Bn-q0B)) + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CB_d*S_CB_d);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / (pi+(1-pi)*rho);
elseif eC == 1
	a_CA_n = 1-ell;
	a_CB_n = 0;
	a_CA_d = 1-ell;
	a_CB_d = 0;
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	S_CB_d = theta*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
	L_A    = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
	L_B_rho_bar = 0;
	S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
	S_CB_tilde  = - i*q0B + ell*(u(q0B,sigma)-q0B);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / (pi+(1-pi)*rho);
end

function [eC,L_A,L_B,q0A,q1A,q0B,q1Bn,q1Bd,eNn,eNd]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% general equilibrium

% initial guess
eC = A/(A+B);

cnt = 0;
precision = 1;
while precision > 10^-6
    cnt = cnt + 1;
    if cnt > 1
        eC = temp_eC;
    end
    
    fprintf('finding a fixed point ... iter: %.0f ... ell = %.2f, theta = %.2f, eta = %.2f, i = %.4f, pi = %.5f, A= %.4f, B = %.4f ... eC = %.4f ',cnt,ell,theta,eta,i,pi,A,B,eC)
    [G_eC,L_A,L_B,q0A,q1A,q0B,q1Bn,q1Bd,eNn,eNd]=partialEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    fprintf('... G_eC = %.4f\n\n',G_eC)
    
    itr = 200; %1
    if cnt < itr
        if eC < 1
            if eC > 0
%                 temp_eC = max(0,min(1,eC + .5*G_eC)); 
                temp_eC = max(0,min(1,eC + 5*G_eC)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        precision = abs(G_eC);
    elseif cnt >= itr
        if cnt > itr
            temp_eC_2 = temp_eC_1;
        end
        if eC < 1
            if eC > 0
                temp_eC_1 = eC;
                temp_eC = max(0,min(1,eC + sign(G_eC)*10^-4)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        if cnt > itr
            if temp_eC_2 == temp_eC
                break
            end
        end
    end
end

%=========================================================================
% calibration
%-------------------------------------------------------------------------
function D=distance(y,ell,theta,sigma,rho,i,M,A_mat,B_mat,TARGET)
% measure the distance of moments from targets

ETA = y(1);
PI  = y(2);

J_A   = zeros(1,length(A_mat));
J_B   = zeros(1,length(A_mat));
for j=1:length(A_mat)
    A = A_mat(j);
    B = B_mat(j);
    [~,L_A,L_B,~,~,~,~,~,~,~]=fixedPoint(ell,theta,sigma,PI,rho,ETA,i,M,A,B);
	pi_bar = PI+(1-PI)*rho;
	J_A(j) = (i*100 - L_A)*4;
	J_B(j) = (i*100 - L_B + (1/pi_bar-1)*100)*4;
end
OUTPUT = [J_A(1)        J_A(2);
          J_B(1)        J_B(2);
          J_B(1)-J_A(1) J_B(2)-J_A(2)];
D = (TARGET(3,1)-OUTPUT(3,1)).^2 + (TARGET(3,2)-OUTPUT(3,2)).^2;

function [ell,theta,sigma,pi,rho,eta,i,M,A_mat,B_mat]=calibration(theta,M)

% (other) calibrated parameters
ell = 0.5; sigma = 0.34; rho = 0.4; i = 0.0175;

% 2004 to 2011; 2012 to 2018
A_mat  = [0.0109 0.0049];
B_mat  = [0.0477 0.0445];
%          pre     post
TARGET = [4.3555  2.6408;  % j_A
          4.5801  2.5999;  % j_B
          0.2245 -0.0409]; % dif

% calibrate pi and eta, for a given theta
fun = @(y) distance(y,ell,theta,sigma,rho,i,M,A_mat,B_mat,TARGET);
lb = [0,   0.99]; % lower bound constraint
ub = [0.3, 1]; % upper bound constraint
% lb = [0,   0.994]; % lower bound constraint
% ub = [0.3, 0.996]; % upper bound constraint
y0 = (lb + ub)/2;  % initial guess
options = optimoptions('fmincon','Display','off','FunctionTolerance',1e-10);
y = fmincon(fun,y0,[],[],[],[],lb,ub,[],options);
eta = y(1);
pi  = y(2);

function [theta,pi,eta,L_A,L_B,J_A,J_B,DIF,DID,Q0A,Q0B]=details(ell,theta,sigma,pi,rho,eta,i,M,A_mat,B_mat)

L_A = zeros(1,length(A_mat));
L_B = zeros(1,length(A_mat));
J_A = zeros(1,length(A_mat));
J_B = zeros(1,length(A_mat));
Q0A = zeros(1,length(A_mat));
Q0B = zeros(1,length(A_mat));
for j=1:length(A_mat) 
    A = A_mat(j);
    B = B_mat(j);
    [~,L_A(j),L_B(j),Q0A(j),~,Q0B(j),~,~,~,~]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B);
    pi_bar = pi+(1-pi)*rho;
    J_A(j) = (i*100 - L_A(j))*4;
    J_B(j) = (i*100 - L_B(j) + (1/pi_bar-1)*100)*4;
end
DIF = J_B-J_A;
DID = DIF(1)-DIF(2);

%=========================================================================
% tables
%-------------------------------------------------------------------------
function table2(ell,theta,sigma,pi,rho,eta,i,M,A,B)

% calibrated parameters (table 1)
pi = 0.996; rho = 0.4; sigma = 0.34; ell = 0.5; theta = 0.8; eta = 0.10; i = 0.0175;

% 2004 to 2011; 2012 to 2018
A_mat  = [0.0109 0.0049];
B_mat  = [0.0477 0.0445];

[~,~,~,~,~,~,~,DIF,DID,Q0A,Q0B]=details(ell,theta,sigma,pi,rho,eta,i,M,A_mat,B_mat);

PrePeriod  = round(DIF(1),4);
PostPeriod = round(DIF(2),4);
Change     = round(DID,4);
table(PrePeriod,PostPeriod,Change)

q0 = ["q0A";"q0B"];
PrePeriod  = round([Q0A(1);Q0B(1)],4);
PostPeriod = round([Q0A(2);Q0B(2)],4);
table(q0,PrePeriod,PostPeriod)

function table3(ell,theta,sigma,pi,rho,eta,i,M,A,B)

Theta = [0.6;0.7;0.8;0.9;0.99];
Pi    = zeros(length(Theta),1);
Eta   = zeros(length(Theta),1);
for j=1:length(Theta)
    [~,~,~,pi,~,eta,~,~,~,~]=calibration(Theta(j),M);
    Pi(j)  = round(pi,4);
    Eta(j) = round(eta,3);
end
table(Theta,Pi,Eta)
